from libmyatomism import *
#from math import *
#import array 
#import matplotlib.pyplot as plt

msLogger.setPriorityMin(Priority.ERROR)

unit=msUnitsManager.New("Angstrom Degree amu kcal/mol ps")

q0 = msScalarVariable.New("Angstrom").set(5,0.5,6,0.1,0.1).setId("q0")
q1 = msScalarVariable.New("Degree").set(120,0,360,2,2).setId("q1")
q2 = msScalarVariable.New("Degree").set(120,0,360,2,2).setId("q2")
q3 = msScalarVariable.New("Degree").set(120,0,360,2,2).setId("q3")

gcoors = msGeneralizedCoordinates.New(unit).addVar(q0).addVar(q1).addVar(q2).addVar(q3)

zmat = msZMat.New(unit).setId("zmat").set("X N 0 0.55 N 0 0.55 1 179.9 X 0 5. 1 90. 2 0. N 3 0.55 0 90. 1 180. N 3 0.55 4 179.9 0 180.")

kinfct0 = msParser.New("Angstrom").setCoordinates(gcoors).setExpression("q0").setId("f=q0")
kinfct1 = msParser.New("Degree").setCoordinates(gcoors).setExpression("q1").setId("f=q1")
kinfct2 = msParser.New("Degree").setCoordinates(gcoors).setExpression("q2").setId("f=q2")
kinfct3 = msParser.New("Degree").setCoordinates(gcoors).setExpression("q3").setId("f=q3")

system=msSystem.New(unit).addEntity(zmat)

kinop=msKineticOperator.New(unit).set(gcoors,system)
kinop.addDynamicDof(zmat,3,kinfct0).addDynamicDof(zmat,4,kinfct1).addDynamicDof(zmat,7,kinfct2).addDynamicDof(zmat,8,kinfct3)

cfour = qcCfour.New().setLevel( BasisSet.ANO0 , CiType.CCSDpT ).setGenBasSymLink("/org/centers/pecos/APPLICATIONS/cfour/basis/GENBAS")

zmatN2 = msZMat.New(unit).setId("zmat").set("N N 0 1.1")
eN2    = cfour.singlePoint(zmatN2)

#scenarios=[[0.1,0.1,0],[90.,0.1,0],[90,90,0],[90,90,90]]
scenarios=[[45,0.1,0],[45,45,0],[45,0.1,45],[0.1,45,45]]
for scenar in scenarios:
    filename="approach_"+str(scenar[0])+"_"+str(scenar[1])+"_"+str(scenar[2])+".txt"
    file = open(filename,'w')
    gcoors.setCoordinate(1,scenar[0]).setCoordinate(2,scenar[1]).setCoordinate(3,scenar[2])

    for i in range(0,50):
        gcoors.setCoordinate(0,6-4.5*i/50.)
        kinop.setDynamicDofs()
        e = cfour.singlePoint(zmat)
        print q0.Value,"\t",q1.Value,"\t",q2.Value,"\t",q3.Value,"\t",e,"\t",unit.convert("hartree",e-2*eN2)
        file.write(str(q0.Value)+" "+str(q1.Value)+" "+str(q2.Value)+" "+str(q3.Value)+" "+str(e)+" "+str(unit.convert("hartree",e-2*eN2))+"\n" )
    file.close()
    